clear; close all; clc;

saveArg = 0;
listFile = dir('*.dat');

rho = [910*ones(40,1);1000*ones(3,1);...
    1000*ones(35,1);1000*ones(12,1)];
drho = [100*ones(40,1);10*ones(3,1);...
    70*ones(35,1);10*ones(12,1)];
mu = [0.1*ones(40,1);0.001*ones(3,1);...
    0.01*ones(47,1)];
Q = [11*ones(25,1);170*ones(3,1);2*ones(9,1);...
    21.3*ones(3,1);299.3;72.68;36.87;...
    38*ones(35,1);3*ones(12,1)]*10^-6/60;
E = [3620*ones(25,1);3520*ones(3,1);1790*ones(9,1);...
    1420*ones(3,1);1840*ones(3,1);...
    3990*ones(18,1);4580*ones(17,1);2400*ones(12,1)];
g = 9.8;
poisson = 0.5;

frames = [1100:400:10700, 2200,2700,3700, ...
    7500:3500:32000,39000, 4100,10100,20700, 600,2100,4200, ...
    250:1000:17250, 2250:1000:18250, 250:1000:11250]';
t = (frames - [200*ones(25,1);230*ones(3,1);3262*ones(9,1);...
    220*ones(3,1);zeros(3,1);...
    -9520*ones(18,1);zeros(17,1);zeros(12,1)])./...
    ([10*ones(25,1);20*ones(3,1);10*ones(9,1);...
    50*ones(3,1);50*ones(3,1);...
    60*ones(18,1);50*ones(17,1);ones(12,1)]);
t(end-11:end) = t(end-11:end)*0.25 + 1013;

frac_data = xlsread('..\fracture_vel.xlsx');
dLdt = frac_data(:,2);
dBdt = frac_data(:,3);
dHdt = frac_data(:,4);
clear frac_data

Vol = Q.*t;
Vol(41:43) = 10^-5*[3.3,14,33.8];
G = E/(2*(1+0.5));
Kc = sqrt(2*1*E);
PMenand = (G(1:43).*(drho(1:43)*g).^2.*Q(1:43)./...
    dLdt/(1-poisson)).^(1/3);

colors = [repmat([0.4,0.7,0],25,1);repmat([0.9,0.9,0],3,1);...
    repmat([0,0,1],9,1);repmat([0.4,0.7,0],3,1);...
    repmat([0.9,0.9,0],3,1);...
    repmat([0.4,0.7,0],35,1);repmat([0,0,1],12,1)];
% colorface = [repmat(1*[1,1,1],43,1);...
%     repmat([0.4,0.7,0],35,1);repmat([0,0,1],12,1)];
colorface = colors;
colorface(44:end,:) = 1;
symbols = [repmat({'o'},37,1);repmat({'d'},6,1);repmat({'o'},47,1)];


%% Read data

x = nan(10^6,length(listFile));
y = x;
u = x;
v = x;
h = waitbar(0,'Loading data...');
for i = 1:length(listFile)
    data = importdata(sprintf('%s',listFile(i).name));
    data = data.data;
    xp = data(:,1);
    yp = data(:,2);
    up = data(:,3);
    vp = data(:,4);
    
    xp = xp(vp~=0);
    yp = yp(vp~=0);
    up = up(vp~=0);
    vp = vp(vp~=0);
    
    x(1:length(xp),i) = xp;
    y(1:length(xp),i) = yp;
    u(1:length(xp),i) = up;
    v(1:length(xp),i) = vp;

    waitbar(i/length(listFile),h);
end
close(h); clear h;

dataSize = sum(u,2,'omitnan');
dataSize = find(dataSize~=0,1,'last');
x = x(1:dataSize,:);
y = y(1:dataSize,:);
u = u(1:dataSize,:);
v = v(1:dataSize,:);
clear data* *p;

% Adjust x & y pos
x = x - [15*ones(size(x,1),25),10*ones(size(x,1),3),...
    5*ones(size(x,1),9),5*ones(size(x,1),3),40*ones(size(x,1),3),...
    0*ones(size(x,1),18),20*ones(size(x,1),17),-20*ones(size(x,1),12)];
y(:,1:25) = y(:,1:25) - min(y(:,1:25),[],'all');
y(:,26:28) = y(:,26:28) - min(y(:,26:28),[],'all');
y(:,29:37) = y(:,29:37) - min(y(:,29:37),[],'all');
y(:,38:40) = y(:,38:40) - min(y(:,38:40),[],'all');
y(:,41:43) = y(:,41:43) - min(y(:,41:43),[],'all');
y(:,44:61) = y(:,44:61) - min(y(:,44:61),[],'all');
y(:,62:78) = y(:,62:78) - min(y(:,62:78),[],'all');
y(:,79:90) = y(:,79:90) - min(y(:,79:90),[],'all');
x = x/1000;
y = y/1000;

L = range(y,1)';
B = range(x,1)';
H = Vol./(pi/4*L.*B);
B(79:90) = max(x(:,79:90))*2;


Pb = drho*g.*L;
PbB = drho*g.*B;
Pf = Kc./sqrt(B);
Pe = G.*H./B/(1-poisson);
Pv = 12*mu.*Q.*L./B./H.^3;
Hef = Kc*(1-poisson)./G.*sqrt(2*B/pi);
Re = rho.*Q.*H./mu./L./B;
% dPMen = ((G/(1-poisson)).*(drho*g).^2.*Q./dLdt*pi).^(1/3);


%% 1st check

check = 90;

figure(1000); clf;
set(gcf,'units','normalized','OuterPosition',[0.1 0.2 0.6 0.6]);
subplot(1,2,1);
    scatter(x(:,check),y(:,check),30,u(:,check),'filled');
    colorbar;
%     caxis(max(abs(u(:,check)),[],'all')*0.1*[-1,1]);
    ylim([0,max(y,[],'all')]);
subplot(1,2,2);
    scatter(x(:,check),y(:,check),30,v(:,check),'filled');
    colorbar;
%     caxis(max(abs(v(:,check)),[],'all')*0.1*[-1,1]);
    ylim([0,max(y,[],'all')]);


%% Filtering

stencil = 3;
filtThresh = 10^1;

xs = floor(min(x,[],'all')*100)/100:0.002:ceil(max(x,[],'all')*100)/100;
ys = floor(min(y,[],'all')*100)/100:0.002:ceil(max(y,[],'all')*100)/100;
uMed = nan(length(ys),length(xs),length(listFile));
vMed = uMed;
uMad = uMed;
vMad = uMed;
uMed1 = uMed;
vMed1 = uMed;
vectorMed = uMed;
vectorMad = uMed;
thetaMed = uMed;
thetaMad = uMed;
h = waitbar(0,'Filtering...');
for i1 = 1:length(listFile)
indStartX = find(xs>=min(x(:,i1)),1,'first');
indEndX = find(xs<=max(x(:,i1)),1,'last');
indStartY = find(ys>=min(y(:,i1)),1,'first');
indEndY = find(ys<=max(y(:,i1)),1,'last');
for i2 = indStartY:indEndY
for i3 = indStartX:indEndX
    indY = i2-1:i2+1;
    indY(indY<indStartY) = [];
    indY(indY>indEndY) = [];
    indX = i3-1:i3+1;
    indX(indX<indStartX) = [];
    indX(indX>indEndX) = [];
    indFilt = x(:,i1)>=xs(indX(1)) & x(:,i1)<=xs(indX(end)) & ...
        y(:,i1)>=ys(indY(1)) & y(:,i1)<=ys(indY(end)) & ...
        isnan(u(:,i1))==0;
        
    if sum(indFilt)>0
        indY = i2-stencil:i2+stencil;
        indY(indY<indStartY) = [];
        indY(indY>indEndY) = [];
        indX = i3-stencil:i3+stencil;
        indX(indX<indStartX) = [];
        indX(indX>indEndX) = [];
        indFilt = x(:,i1)>=xs(indX(1)) & x(:,i1)<=xs(indX(end)) & ...
            y(:,i1)>=ys(indY(1)) & y(:,i1)<=ys(indY(end)) & ...
            isnan(u(:,i1))==0;
        
        
        subset = u(indFilt,i1);
        uMed1(i2,i3,i1) = median(subset);
        uMad(i2,i3,i1) = mad(subset);
        subset = v(indFilt,i1);
        vMed1(i2,i3,i1) = median(subset);
        vMad(i2,i3,i1) = mad(subset);
        subset = acos(v(indFilt,i1)./...
            (u(indFilt,i1).^2+v(indFilt,i1).^2).^(1/2));
        thetaMed(i2,i3,i1) = median(subset);
        thetaMad(i2,i3,i1) = mad(subset);
        subset = (u(indFilt,i1).^2+v(indFilt,i1).^2).^(1/2);
        vectorMed(i2,i3,i1) = median(subset);
        vectorMad(i2,i3,i1) = mad(subset);
        
        if abs(uMad(i2,i3,i1)./uMed1(i2,i3,i1)) * ...
                abs(vMad(i2,i3,i1)./vMed1(i2,i3,i1)) < filtThresh
            
            subset = u(indFilt,i1);
            submed = median(subset);
            submad = mad(subset);
            subset(subset<submed-submad|subset>submed+submad) = [];
            uMed(i2,i3,i1) = median(subset);

            subset = v(indFilt,i1);
            submed = median(subset);
            submad = mad(subset);
            subset(subset<submed-submad|subset>submed+submad) = [];
            vMed(i2,i3,i1) = median(subset);
        end
    end
end
waitbar((i1-1)/length(listFile) + ...
(i2-1)/length(ys)/length(listFile),h);
end
end
close(h); clear h;

x = xs;
y = ys;
u = uMed;
v = vMed;

clear i* sub* gridSize xs ys xc yc uMed vMed;


%% 2nd check

Lb = (Kc./drho/g).^(2/3);

% figure(1001); clf;
% for i = 1:length(listFile)
%     pcolor(log10(abs(uMad(:,:,i)./uMed1(:,:,i)).*...
%         abs(vMad(:,:,i)./vMed1(:,:,i))));
%     shading flat; colorbar;
%     caxis([-1,log10(filtThresh)]);
% %     caxis([-1,0.5]);
%     title(sprintf('%d',i));
%     pause();
% end

figure(1002); clf;
for i = 1:length(listFile)
    pcolor(x,y,v(:,:,i)); shading flat; colorbar;
    caxis(max(v(:,:,i),[],'all')*[-1,1]);
    colormap jet;
    title(sprintf('%0.2f',L(i)/Lb(i)));
    pause();
end


%% Vorticity

stencil = 1;

vort = nan(size(u));
h = waitbar(0,'Finding vorticity...');
for i = 1:size(u,3)
    indEndY = sum(u(:,:,i),2,'omitnan');
    indStartY = find(indEndY~=0,1,'first');
    indEndY = find(indEndY~=0,1,'last');
    indEndX = sum(u(:,:,i),1,'omitnan');
    indStartX = find(indEndX~=0,1,'first');
    indEndX = find(indEndX~=0,1,'last');
for j1 = indStartY:indEndY
for j2 = indStartX:indEndX
    indY = j1-1:j1+1;
    indY(indY<indStartY) = [];
    indY(indY>indEndY) = [];
    indY(isnan(u(indY,j2,i))) = [];
    indX = j2-1:j2+1;
    indX(indX<indStartX) = [];
    indX(indX>indEndX) = [];
    indX(isnan(v(j1,indX,i))) = [];
    
    if length(indY)>1 && length(indX)>1
        
    indY = j1-stencil:j1+stencil;
    indY(indY<indStartY) = [];
    indY(indY>indEndY) = [];
    indY(isnan(u(indY,j2,i))) = [];
    indX = j2-stencil:j2+stencil;
    indX(indX<indStartX) = [];
    indX(indX>indEndX) = [];
    indX(isnan(v(j1,indX,i))) = [];
    
    if length(indY)>1 && length(indX)>1
        P = polyfit(y(indY)',u(indY,j2,i),1);
        dudy = P(1);
        P = polyfit(x(indX),v(j1,indX,i),1);
        dvdx = P(1);
        vort(j1,j2,i) = dvdx - dudy;
    end
    end
end
end
waitbar((i-1)/(size(u,3)-1),h);
end
close(h); clear h;

clear i* j* stencil dudy dvdx P;

% figure(2000);
% pcolor(vort(:,:,41)); shading flat; colorbar;


%% Cut out left side S3

u(:,1:floor(size(u,2)/2),79:end) = NaN;
v(:,1:floor(size(u,2)/2),79:end) = NaN;
vort(:,1:floor(size(u,2)/2),79:end) = NaN;


%% Cut out right side C1

u(:,76:end,12:25) = NaN;
v(:,76:end,12:25) = NaN;
vort(:,76:end,12:25) = NaN;


%% Pixel flux

dx = x(2)-x(1);
exchParam = nan(length(listFile),1);
e2 = exchParam;
for i = 1:length(listFile)
    check = v(:,:,i);
    e2(i) = -sum(check(check<0),'all','omitnan')./...
        sum(check(check>0),'all','omitnan');
    if i <= 25
%         maxInd = find(check==max(check,[],'all'));
%         [maxY,maxX] = ind2sub(size(check),maxInd);
        topInd = sum(check,2,'omitnan');
        topInd = find(topInd~=0,1,'last');
        leftInd = find(~isnan(check(topInd,:)),1,'first');
        rightInd = find(~isnan(check(topInd,:)),1,'last');
        centInd = round(mean([leftInd,rightInd]),0);
        
        check(:,centInd+1:end) = NaN;
    end

%     topInd = sum(check,2,'omitnan');
%     topInd = find(topInd~=0,1,'last');
%     leftInd = find(~isnan(check(topInd,:)),1,'first');
%     rightInd = find(~isnan(check(topInd,:)),1,'last');
%     gridY = y(1:topInd);
%     gridY = gridY-mean(gridY);
%     gridY = sqrt(1-(gridY/max(gridY)).^2).*(1+gridY/max(gridY));
%     gridX = x(leftInd:rightInd);
%     gridX = sqrt(1-(gridX/max(abs(gridX))).^2);
%     gridH = real(gridY'*gridX);
%     gridVol = sum(gridH,'all')*(x(2)-x(1))*(y(2)-y(1));
%     gridH = gridH/gridVol*Vol(i);
    
    fluxUp = check;
    fluxUp(fluxUp<0) = 0;
    fluxUp = sum(fluxUp*H(i)*dx,2,'omitnan');
    fluxDown = check;
    fluxDown(fluxDown>0) = 0;
    fluxDown = sum(fluxDown*H(i)*dx,2,'omitnan');
    
    zeroInd = fluxUp==0 | fluxDown==0;
    fluxUp(zeroInd) = [];
    fluxDown(zeroInd) = [];
    
    exchParam(i) = abs(mean(fluxDown))./mean(fluxUp);
%     exchParam(i) = e2(i);
    
%     figure(1000);
%     pcolor(check); shading flat; colorbar;
%     pause();
end

xr = PbB./Pf;
yr = exchParam;
figure(10001); clf;
for i = 1:length(xr)
    plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',colors(i,:),...
        'MarkerFaceColor',colorface(i,:)); hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
grid on;
xlabel('\Pi_P = P_b/P_f');
ylabel('\Pi_v');

yr(44:end) = NaN;
modelA = logspace(0.5,1.5,30);
modelB = logspace(0,1,30);
modelC = logspace(-3,-2,30);
errorMat = nan(length(modelA),length(modelB),length(modelC));
for m1 = 1:length(modelA)
for m2 = 1:length(modelB)
for m3 = 1:length(modelC)
    model = testeq(modelA(m1),modelB(m2),modelC(m3),xr);
    SSR = sum((log10(model) - log10(yr)).^2,'omitnan');
    SST = sum((log10(model) - log10(mean(yr,'omitnan'))).^2,'omitnan');
    R2 = 1-SSR/SST;
    errorMat(m1,m2,m3) = R2;
end
end
end
minInd = find(errorMat==max(errorMat,[],'all'));
[indA,indB,indC] = ind2sub(size(errorMat),minInd);
bestA = modelA(indA);
bestB = modelB(indB);
bestC = modelC(indC);
model = testeq(bestA,bestB,bestC,xr);

surfx = logspace(log10(min(xr)),log10(max(xr)),100);
surfy = testeq(bestA,bestB,bestC,surfx);
plot(surfx,surfy,'k--');
% ylim([10^-3,1]);
tx = text(0.1,0.05,...
    sprintf('\\Pi_v = exp(-%0.1f\\Pi_P^{%0.1f})+%0.4f',...
    bestA,bestB,bestC));
set(tx,'Fontsize',14);
tx = text(0.1,0.02,...
    sprintf('R^2 = %0.2f',errorMat(indA,indB,indC)));
set(tx,'Fontsize',14);

% plot(xr,e2,'rx');


%% Cumulative density function fitting

Lb = (Kc./drho/g).^(2/3);

cumMed = nan(length(listFile),1);
cumMad = cumMed;
figure(19); clf;
for i = 1:90%length(listFile)
    check = v(:,:,i);
    check = check(~isnan(check));
    bounds = minmax(check,[0.05 0.95]);
    check(check<bounds(1)|check>bounds(2)) = [];
    xs = linspace(min(check),max(check),300);
    
    cumMed(i) = median(check,'omitnan');
    cumMad(i) = mad(check,1);
    fnorm = cdf('Normal',xs,cumMed(i),cumMad(i));
    
%     figure(18); clf;
%     cp = cdfplot(check); hold on;
%     cp.Color = 'k';
%     plot(xs,fnorm,'k--');
%     pause();
    
    figure(19);
    plot(L(i)./Lb(i),cumMad(i)./cumMed(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
% set(gca,'XScale','log')
% set(gca,'YScale','log');


%% Dimensionless exchange

exchMed = nan(length(listFile),1);
exchMad = exchMed;
medVort = exchMed;
medVel = exchMed;
medU = exchMed;
medVec = exchMed;
meanVort = exchMed;
meanVel = exchMed;
meanU = exchMed;
meanVec = exchMed;
madVort = exchMed;
madVel = exchMed;
madU = exchMed;
madVec = exchMed;
devVort = exchMed;
devVel = exchMed;
devU = exchMed;
devVec = exchMed;
maxVel = exchMed;
maxVec = exchMed;
minVel = exchMed;
minVec = exchMed;

dimExch = abs(vort)./v;
vector = (u.^2+v.^2).^(1/2);

for i = 1:length(listFile)
    check = v(:,:,i);
%     if i>=12 && i<=25
%         check(:,76:end) = NaN;
%     end
    medVel(i) = median(check,'all','omitnan');
    meanVel(i) = mean(check,'all','omitnan');
    madVel(i) = mad(check,1,'all');
    devVel(i) = std(check,0,'all','omitnan');
    maxVel(i) = max(check,[],'all');
    minVel(i) = min(check,[],'all');
    
    check = u(:,:,i);
%     if i>=12 && i<=25
%         check(:,76:end) = NaN;
%     end
    medU(i) = median(abs(check),'all','omitnan');
    meanU(i) = mean(abs(check),'all','omitnan');
    madU(i) = mad(abs(check),1,'all');
    devU(i) = std(abs(check),0,'all','omitnan');
    
    check = vector(:,:,i);
%     if i>=12 && i<=25
%         check(:,76:end) = NaN;
%     end
    medVec(i) = median(check,'all','omitnan');
    meanVec(i) = mean(check,'all','omitnan');
    madVec(i) = mad(check,1,'all');
    devVec(i) = std(check,0,'all','omitnan');
    maxVec(i) = max(check,[],'all');
    minVec(i) = min(check,[],'all');
    
    check = vort(:,:,i);
%     if i>=12 && i<=25
%         check(:,76:end) = NaN;
%     end
    medVort(i) = median(abs(check),'all','omitnan');
    meanVort(i) = mean(abs(check),'all','omitnan');
    madVort(i) = mad(abs(check),1,'all');
    devVort(i) = std(abs(check),0,'all','omitnan');
end
clear exchS* i;


%% Plot dimensionless exchange

figure(10); clf;
set(gcf,'Units','normalized','outerposition',[0 0.05 1 0.95]);
k = 0;
for i = [1,26,29,38,41, 10,27,33,39,42, 25,28,37,40,43]
    k = k + 1;
    subplot(3,5,k);
    pcolor(x,y,abs(dimExch(:,:,i)));
    set(gca,'ColorScale','log');
%     caxis([min(abs(dimExch),[],'all'),max(abs(dimExch),[],'all')]);
end
for i = 1:15
    subplot(3,5,i);
    shading flat;
    colorbar;
%     if i == 1 || i == 6 || i == 11
%         caxis([-10^2,10^2]);
%     elseif i == 2 || i == 7 || i == 12
%         caxis([-10^2,10^2]);
%     elseif i == 3 || i == 8 || i == 13
%         caxis([-10^2,10^2]);
%     elseif i == 4 || i == 9 || i == 14
%         caxis([-10^2,10^2]);
%     else
%         caxis([-10^2,10^2]);
%     end
%     caxis([10^-3,10^3]);
%     caxis(10^2*[-1,1]);
    caxis([10^0,10^3]);
    ylim([0,max(y)]);
    xlim([min(x),max(x)]);
end
clear i k;
subplot(3,5,1); ylabel('initial stage'); title('Mid flux, oil');
subplot(3,5,2); title('High flux, oil');
subplot(3,5,3); title('Low flux, oil');
subplot(3,5,4); title('Mid flux, oil');
subplot(3,5,5); title('High flux, water');
subplot(3,5,6); ylabel('middle stage');
subplot(3,5,11); ylabel('late stage');


%% Plot vorticity fields

figure(11); clf;
set(gcf,'Units','normalized','outerposition',[0 0.05 1 0.95]);
k = 0;
for i = [1,26,29,38,41, 10,27,33,39,42, 25,28,37,40,43]
    k = k + 1;
    subplot(3,5,k);
    pcolor(x,y,abs(vort(:,:,i)));
    set(gca,'ColorScale','log');
end
for i = 1:15
    subplot(3,5,i);
    shading flat;
    colorbar;
    if i == 1 || i == 6 || i == 11
        caxis([10^-4,10^-1]);
    elseif i == 2 || i == 7 || i == 12
        caxis([10^-3,10^0]);
    elseif i == 3 || i == 8 || i == 13
        caxis([10^-5,10^-2]);
    elseif i == 4 || i == 9 || i == 14
        caxis([10^-3,10^0]);
    else
        caxis([10^-3,10^1]);
    end
%     caxis([10^-5,10^0]);
    ylim([0,max(y,[],'all')]);
    xlim([min(x,[],'all'),max(x,[],'all')]);
end
clear i k;


%% Plot velocity fields

figure(12); clf;
set(gcf,'Units','normalized','outerposition',[0 0.05 1 0.95]);
k = 0;
for i = [1,26,29,38,41, 10,27,33,39,42, 25,28,37,40,43]
    k = k + 1;
    subplot(3,5,k);
    pcolor(x,y,vector(:,:,i));
    set(gca,'ColorScale','log');
end
for i = 1:15
    subplot(3,5,i);
    shading flat;
    colorbar;
    if i == 1 || i == 6 || i == 11
        caxis([10^-4,10^-2]);
    elseif i == 2 || i == 7 || i == 12
        caxis([10^-4,10^-2]);
    elseif i == 3 || i == 8 || i == 13
        caxis([10^-6,10^-4]);
    elseif i == 4 || i == 9 || i == 14
        caxis([10^-4,10^-2]);
    else
        caxis([10^-4,10^-1]);
    end
%     caxis(10^-2*[-1,1]);
    ylim([0,max(y,[],'all')]);
    xlim([min(x,[],'all'),max(x,[],'all')]);
end
clear i k;


%% Plot stats pixel flux

xr = PbB./Pf;
yr = exchParam;

figure(20); clf;
set(gcf,'Units','normalized','OuterPosition',[0.1 0.2 0.8 0.7]);
subplot(1,2,1);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'kd','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','k','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0,0,1],...
    'MarkerEdgeColor',[0,0,1],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
for i = 1:length(xr)
    if i <= 43
        cval = colors(i,:);
        cvalf = cval;
    elseif i <= 78
        cval = [0.85 0.925 0.75];
        cvalf = 'none';
    else
        cval = [0.8 0.8 1];
        cvalf = 'none';
    end
    pl = plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',cval,'markerfacecolor',cvalf,...
        'LineWidth',2);
    hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
set(gca,'FontSize',16);
grid on;
xlabel('\Pi_P = P_b/P_f');
ylabel('\Pi_Q = |Q_d|/Q_a');
xlimits = [0.02,2];
ylimits = [10^-3,5];
xlim(xlimits);
ylim(ylimits);
tx = text(10^(log10(xlimits(1))+0.92*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'a');
set(tx,'FontSize',24);
plot(get(gca,'XLim'),[1,1],'k--');
plot(get(gca,'XLim'),[1,1],'k--');
tx = text(0.025,0.65,'circulating pattern');
set(tx,'Fontsize',16);
% tx = text(0.09,0.0017,'ascending pattern \rightarrow');
tx = text(0.15,0.0017,'ascending pattern');
set(tx,'Fontsize',16);

% leg = legend('pump','gravity','C series','S series',...
%     'low flux','med flux','high flux',...
%     'NumColumns',2,'fontsize',12,'location','southwest');
% set(leg,'Position',[0.11 0.30 0.22 0.21]);

% yr(44:end) = NaN;
% modelA = logspace(0.5,1.5,30);
% modelB = logspace(0,1,30);
% modelC = logspace(-3,-2,30);
% errorMat = nan(length(modelA),length(modelB),length(modelC));
% for m1 = 1:length(modelA)
% for m2 = 1:length(modelB)
% for m3 = 1:length(modelC)
%     model = testeq(modelA(m1),modelB(m2),modelC(m3),xr);
%     SSR = sum((log10(model) - log10(yr)).^2,'omitnan');
%     SST = sum((log10(model) - log10(mean(yr,'omitnan'))).^2,'omitnan');
%     R2 = 1-SSR/SST;
%     errorMat(m1,m2,m3) = R2;
% end
% end
% end
% minInd = find(errorMat==max(errorMat,[],'all'));
% [indA,indB,indC] = ind2sub(size(errorMat),minInd);
% bestA = modelA(indA);
% bestB = modelB(indB);
% bestC = modelC(indC);
% surfx = logspace(log10(xlimits(1)),log10(xlimits(2)),100);
% surfy = testeq(bestA,bestB,bestC,surfx);
% plot(surfx,surfy,'k--');
% tx = text(0.035,0.03,...
%     sprintf('\\Pi_Q = exp(-%0.1f\\Pi_P^{%0.1f})+%0.4f',...
%     bestA,bestB,bestC));
% set(tx,'Fontsize',14);
% tx = text(0.035,0.01,...
%     sprintf('R^2 = %0.2f',errorMat(indA,indB,indC)));
% set(tx,'Fontsize',14);

xr = Re;
yr = exchParam;
subplot(1,2,2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'kd','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','k','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0,0,1],...
    'MarkerEdgeColor',[0,0,1],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
for i = 1:length(xr)
    if i <= 43
        cval = colors(i,:);
        cvalf = cval;
    elseif i <= 78
        cval = [0.85 0.925 0.75];
        cvalf = 'none';
    else
        cval = [0.8 0.8 1];
        cvalf = 'none';
    end
    pl = plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',cval,'markerfacecolor',cvalf,...
        'LineWidth',2);
    hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
xlimits = [6e-5,10];
ylimits = [10^-3,5];
xlim(xlimits);
ylim(ylimits);
tx = text(10^(log10(xlimits(1))+0.92*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'b');
set(tx,'FontSize',24);
set(gca,'FontSize',16);
grid on;
xlabel('Re');
set(gca,'YTickLabel',[]);
set(gca,'XTick',10.^(-4:2:0));

leg = legend('pump','gravity','C series','S series',...
    'low flux','med flux','high flux',...
    'NumColumns',2,'fontsize',12,'location','southeast');
% set(leg,'Position',[0.11 0.30 0.22 0.21]);

x1 = 0.10; x2 = 0.55; y1 = 0.21;
wd = 0.42; ht = 0.73;
subplot(1,2,1); set(gca,'Position',[x1 y1 wd ht]);
subplot(1,2,2); set(gca,'Position',[x2 y1 wd ht]);


%% Plot model velocity

xr = PbB./Pf; xlab = '\Pi_P = P_b/P_f';
% xr = Pb./Pe; xlab = '\Pi_P = P_b/P_e';
% xr = Pv./Pf; xlab = 'P_v/P_f';
% xr = dP./Pb; xlab = '\DeltaP/P_b';
% xr = Re; xlab = 'Re';
yr = medVec;
% yr2 = Q./H./(B*pi);
% vInflux = Q.*min([L,B],[],2)./Vol;
vInflux = Q./H./(2*pi*sqrt(((B/2).^2+(L/2).^2)/2));
vLister = ((Q./B).^2.*drho.*g./mu/12).^(1/3);
yr4 = Q./pi*4/0.02^2;

figure(30); clf;
set(gcf,'Units','normalized','OuterPosition',[0 0.05 0.4 0.95]);
subplot(5,1,1);
plot(NaN,NaN,'ko','MarkerSize',8,...
    'markerfacecolor','k','LineWidth',2,...
    'LineStyle','none'); hold on;
plot(NaN,NaN,'kd','MarkerSize',8,...
    'markerfacecolor','k','LineWidth',2,...
    'LineStyle','none');
plot(NaN,NaN,'k^','MarkerSize',8,...
    'markerfacecolor','k','LineWidth',2,...
    'LineStyle','none');
plot(NaN,NaN,'kv','MarkerSize',8,...
    'markerfacecolor','k','LineWidth',2,...
    'LineStyle','none');
for i = 1:43
    if i <= 25
        k = 1;
    elseif i <= 28
        k = 2;
    elseif i <= 37
        k = 3;
    elseif i <= 40
        k = 4;
    elseif i <= 43
        k = 5;
    end
    subplot(5,1,k);
    plot(xr(i),yr(i),'marker',symbols{i},'MarkerSize',8,...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
    plot(xr(i),vInflux(i),'v','MarkerSize',4,...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2);
    plot(xr(i),vLister(i),'^','MarkerSize',4,...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2);
    
    set(gca,'XScale','log');
    set(gca,'YScale','log');
    xlim([min(xr),max(xr)]);
end
% ylim([-0.05,1]);

for i = 1:5
    subplot(5,1,i);
    grid on;
    set(gca,'Fontsize',14);
    if i <= 4
        set(gca,'XTickLabel',[]);
    end
    set(gca,'YTick',10.^(-5:-2));
end

subplot(5,1,3);
ylabel('$\bar{\bf{v}}_{med}$ (m/s)','Interpreter','latex');
subplot(5,1,5); xlabel(xlab);
subplot(5,1,1);
v = [0.2, 3e-4; 0.3,2.5e-4; ...
    0.279 2.3e-4; 0.28 2.9e-4];
plot(v(1:2,1),v(1:2,2),'k-','LineWidth',1.5);
plot(v(2:3,1),v(2:3,2),'k-','LineWidth',1.5);
plot(v([2,4],1),v([2,4],2),'k-','LineWidth',1.5);
v = [0.42, 2.7e-4; 0.55,5e-4; ...
    0.51 4.7e-4; 0.52 3.8e-4];
plot(v(1:2,1),v(1:2,2),'k-','LineWidth',1.5);
plot(v(2:3,1),v(2:3,2),'k-','LineWidth',1.5);
plot(v([2,4],1),v([2,4],2),'k-','LineWidth',1.5);
tx = text(0.33,3e-4,'transition');
set(tx,'FontSize',12,'Rotation',20);
leg = legend('pump','gravity',...
    'v_{lub}','v_{inf}',...
    'fontsize',10,'location','west');
% xlim([0.1 0.6]);
% ylim([1e-4 6e-4]);

x1 = 0.17; wd = 0.77; ht = 0.13;
y1 = 0.83; y2 = 0.66; y3 = 0.49;
y4 = 0.32; y5 = 0.15;
subplot(5,1,1); set(gca,'Position',[x1 y1 wd ht]);
subplot(5,1,2); set(gca,'Position',[x1 y2 wd ht]);
subplot(5,1,3); set(gca,'Position',[x1 y3 wd ht]);
subplot(5,1,4); set(gca,'Position',[x1 y4 wd ht]);
subplot(5,1,5); set(gca,'Position',[x1 y5 wd ht]);
set(leg,'Position',[0.18 0.82 0.214 0.146]);


%% More model velocity

% xr = Pb./Pf;
% yr = Re;
% % xr = Pv./Pf;
% % yr = dBdt./dLdt;
% zr = medVec;
% zr = (zr-vInflux)./(vLister-vInflux);
% 
% zr(44:end) = NaN;
% 
% modelA = logspace(0,1,10);
% modelB = logspace(-0.5,0.5,10);
% modelC = logspace(-2.5,-1.5,10);
% modelD = logspace(-1,0,10);
% errorMat = nan(length(modelA),length(modelB),length(modelC),...
%     length(modelD));
% for ma = 1:length(modelA)
% for mb = 1:length(modelB)
% for mc = 1:length(modelC)
% for md = 1:length(modelD)
%         model = (exp(xr/modelA(ma))-modelB(mb)).*...
%             (modelC(mc)*log10(yr)+modelD(md));
%         model(model<0) = 10;
%         errorMat(ma,mb,mc,md) = median((log10(model) - log10(zr)).^2,'omitnan');
% end
% end
% end
% end
% minInd = find(errorMat==min(errorMat,[],'all'));
% [indA,indB,indC,indD] = ind2sub(size(errorMat),minInd);
% bestA = modelA(indA);
% bestB = modelB(indB);
% bestC = modelC(indC);
% bestD = modelD(indD);
% 
% 
% figure(31); clf;
% for i = 1:43
%     plot3(xr(i),yr(i),zr(i),'Marker',symbols{i},...
%         'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
%         'LineWidth',2); hold on;
%     
%     set(gca,'XScale','log');
%     set(gca,'YScale','log');
% %     set(gca,'ZScale','log');
% end
% % view(90,0);
% % xlim([min(xr),max(xr)]);
% % xlim([1.8,2.2]);
% zlim([0,1]);
% xlabel('x');
% ylabel('y');
% zlabel('z');
% 
% xs = logspace(log10(min(xr)),log10(max(xr)),10);
% ys = logspace(log10(min(yr)),log10(max(yr)),10);
% [XS,YS] = meshgrid(xs,ys);
% ZS = (exp(XS/bestA)-bestB).*(bestC*log10(YS)+bestD);
% ZS = exp(XS/10)-1;
% s=surf(XS,YS,ZS);
% alpha(s,0.1);
% text(0.1,1,0.8,sprintf('z=[exp(x/%0.1f)-%0.1f][%0.3flog_{10}(y)+%0.1f]',...
%     bestA,bestB,bestC,bestD));
% view(13,15);


%% Even more model velocity

% xr = Pb./Pe;
% yr = Re;
% zr = vInflux./medVec;
% zr(44:end) = NaN;
% 
% figure(32); clf;
% for i = 1:43
%     plot3(xr(i),yr(i),zr(i),'Marker',symbols{i},...
%         'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
%         'LineWidth',2); hold on;
% end
% set(gca,'XScale','log');
% set(gca,'YScale','log');
% set(gca,'ZScale','log');
% grid on;
% % plot(get(gca,'XLim'),get(gca,'XLim'),'k--');
% view(55,9);
% xlabel('P_b/P_e');
% ylabel('Re');
% zlabel('$v_i/\bf{\bar{v}}$','Interpreter','latex');
% 
% 
% modelA = logspace(0,1,10);
% modelB = logspace(0,1,10);
% modelC = logspace(-1,0,10);
% modelD = logspace(-0.5,0.5,10);
% errorMat = nan(length(modelA),length(modelB),length(modelC),...
%     length(modelD));
% for m1 = 1:length(modelA)
% for m2 = 1:length(modelB)
% for m3 = 1:length(modelC)
% for m4 = 1:length(modelD)
%     model = testeqvi(modelA(m1),modelB(m2),modelC(m3),...
%         modelD(m4),xr,yr);
%     SSR = sum((log10(model) - log10(zr)).^2,'omitnan');
%     SST = sum((log10(model) - log10(mean(zr,'omitnan'))).^2,'omitnan');
%     R2 = 1-SSR/SST;
%     errorMat(m1,m2,m3,m4) = R2;
% end
% end
% end
% end
% minInd = find(errorMat==max(errorMat,[],'all'));
% [indA,indB,indC,indD] = ind2sub(size(errorMat),minInd);
% bestA = modelA(indA);
% bestB = modelB(indB);
% bestC = modelC(indC);
% bestD = modelD(indD);
% bestR2 = errorMat(indA,indB,indC,indD);
% model = testeqvi(bestA,bestB,bestC,bestD,xr,yr);
% 
% for i = 1:length(model)
%     plot3(xr(i)*[1,1],yr(i)*[1,1],[zr(i),model(i)],...
%         'k-','linewidth',2);
% end
% 
% surfx = logspace(log10(min(xr)),log10(max(xr)),10);
% surfy = logspace(log10(min(yr)),log10(max(yr)),10);
% [surfx,surfy] = meshgrid(surfx,surfy);
% surfz = testeqvi(bestA,bestB,bestC,bestD,surfx,surfy);
% s=surf(surfx,surfy,surfz);
% alpha(s,0.2);


%% Model U V velocity

% xr = Pb./Pe;
xr = PbB./Pf;
yr = meanU./meanVec;
figure(40); clf;
set(gcf,'Units','normalized','OuterPosition',[0.1 0.1 0.35 0.85]);
subplot(2,1,1);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'kd','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','k','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0,0,1],...
    'MarkerEdgeColor',[0,0,1],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
for i = 1:length(xr)
    plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
set(gca,'Fontsize',14);
set(gca,'XScale','log');
set(gca,'YScale','log');
grid on;
xlimits = [0.03 2];
ylimits = [0.03 1];
xlim(xlimits); ylim(ylimits);
set(gca,'XTickLabel',[]);
ylabel('$\overline{|u|}/\bar{\bf{v}}$','Interpreter','latex');
tx = text(10^(log10(xlimits(1))+0.91*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'a');
set(tx,'Fontsize',20);
% set(gca,'Position',[0.15 0.19 0.80 0.76]);
leg = legend('pump','gravity','C series','S series',...
    'low flux','med flux','high flux',...
    'NumColumns',2,'fontsize',10,'location','southwest');
set(leg,'Position',[0.20 0.61 0.47 0.14]);

% xr = Pb./Pe;
yr = meanVel./meanVec;
subplot(2,1,2);
for i = [1:42,44:length(xr)]
    plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
set(gca,'Fontsize',14);
set(gca,'XScale','log');
% set(gca,'YScale','log');
grid on;
xlimits = [0.03 2];
ylimits = [-0.5 1.05];
xlim(xlimits); ylim(ylimits);
tx = text(10^(log10(xlimits(1))+0.05*range(log10(xlimits))),...
    0.9,'b');
set(tx,'Fontsize',20);
xlabel('\Pi_P = P_b/P_f');
ylabel('$\bar{v}/\bar{\bf{v}}$','Interpreter','latex');

x1 = 0.19; wd = 0.76; ht = 0.36;
y1 = 0.60; y2 = 0.15;
subplot(2,1,1); set(gca,'Position',[x1 y1 wd ht]);
subplot(2,1,2); set(gca,'Position',[x1 y2 wd ht]);


%% Model Vorticity

% xr1 = maxVec*2./B;
xr1 = (maxVel-minVel)./(B/2);
% xr2 = Pb./Pe;
xr2 = PbB./Pf;
% yr = medVort;
yr = meanVort;


figure(50); clf;
set(gcf,'Units','normalized','OuterPosition',[0.1 0.2 0.9 0.6]);
subplot(1,3,1);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'kd','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','k','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'ko','markersize',8,...
    'MarkerFaceColor','w','LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0,0,1],...
    'MarkerEdgeColor',[0,0,1],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
for i = 1:length(xr1)
    plot(xr1(i),yr(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
set(gca,'Fontsize',14);
grid on;
set(gca,'XTick',10.^(-3:1));
xlimits = [0.0004,4];
ylimits = [0.0004,4];
xlim(xlimits);
ylim(ylimits);
plot(get(gca,'XLim'),get(gca,'XLim'),'k--');
xlabel('$2(v_{max}-v_{min})/B (1/s)$','Interpreter','latex');
ylabel('$\overline{|\omega|} (1/s)$','Interpreter','latex');
tx = text(0.3,0.6,'1:1');
set(tx,'FontSize',14,'Rotation',40);
tx = text(10^(log10(xlimits(1))+0.05*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'a');
set(tx,'FontSize',20);
% set(gca,'Position',[0.15 0.19 0.80 0.76]);

leg = legend('pump','gravity','C series','S series',...
    'low flux','med flux','high flux',...
    'NumColumns',2,'fontsize',12,'location','southeast');
set(leg,'Position',[0.17 0.235 0.18 0.22]);


subplot(1,3,2);
for i = 1:length(xr2)
    plot(xr2(i),yr(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
set(gca,'Fontsize',14);
set(gca,'XTick',10.^(-3:1));
grid on;
xlimits = [0.03,2];
ylimits = [0.0004,3];
xlim(xlimits);
ylim(ylimits);
xlabel('\Pi_P = P_b/P_f');
ylabel('$\overline{|\omega|} (1/s)$','Interpreter','latex');
tx = text(10^(log10(xlimits(1))+0.05*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'b');
set(tx,'FontSize',20);
% set(gca,'Position',[0.15 0.19 0.80 0.76]);


subplot(1,3,3);
for i = 1:length(xr2)
    plot(xr2(i),yr(i)./xr1(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
set(gca,'Fontsize',14);
grid on;
set(gca,'XTick',10.^(-3:1));
set(gca,'YTick',[0.2:0.4:1.4]);
xlimits = [0.03,2];
ylimits = [0.2,1.4];
xlim(xlimits);
ylim(ylimits);
xlabel('\Pi_P');
ylabel('\Pi_\omega');
tx = text(10^(log10(xlimits(1))+0.88*range(log10(xlimits))),...
    10^(log10(ylimits(1))+0.92*range(log10(ylimits))),'c');
set(tx,'FontSize',20);

x1 = 0.08; x2 = 0.43; x3 = 0.74;
y1 = 0.23; wd = 0.23; ht = 0.70;
subplot(1,3,1); set(gca,'Position',[x1 y1 wd ht]);
subplot(1,3,2); set(gca,'Position',[x2 y1 wd ht]);
subplot(1,3,3); set(gca,'Position',[x3 y1 wd ht]);


%% Probability density plots

Lb = (Kc./drho/g).^(2/3);
x1 = 0.07;
y1 = 0.60; y2 = 0.14; wd = 0.16; ht = 0.36; sp = 0.027;
x2 = x1+1*(wd+sp); x3 = x1+2*(wd+sp);
x4 = x1+3*(wd+sp); x5 = x1+4*(wd+sp);
xp = [x1 x2 x3 x4 x5 x1 x2 x3];
yp = [y1 y1 y1 y1 y1 y2 y2 y2];
colies = [0;0.6;0.85]*ones(1,3);

figure(61); clf;
set(gcf,'Units','normalized','OuterPosition',[0.1 0.1 0.8 0.8]);
for i = 1:8
    if i == 1
        checkers = [1,13,25];
        titlab = 'a) C1';
    elseif i == 2
        checkers = 26:28;
        titlab = 'b) C2';
    elseif i == 3
        checkers = [29,34,37];
        titlab = 'c) C3';
    elseif i == 4
        checkers = 38:40;
        titlab = 'd) C4';
    elseif i == 5
        checkers = 41:43;
        titlab = 'e) C5';
    elseif i == 6
        checkers = [44,52,61];
        titlab = 'f) S1';
    elseif i == 7
        checkers = [62,69,78];
        titlab = 'g) S2';
    elseif i == 8
        checkers = [79,83,90];
        titlab = 'h) S3';
    end
    subplot(2,5,i);
    numies = 20;
    for k = 1:3
        check = v(:,:,checkers(k))*1000;
        check = check(~isnan(check));
        bounds = minmax(check,[0.05 0.95]);
        check(check<bounds(1)|check>bounds(2)) = [];
        [bnnum,bnloc] = hist(check,numies);
        plot(bnloc,bnnum./max(bnnum),'-',...
            'Color',colies(k,:),...
            'LineWidth',2); hold on;
    end
    xlimies = get(gca,'XLim');
    tx = text(xlimies(1)+range(xlimies)*0.03,1.35,titlab);
    set(tx,'FontSize',14);
    lg = legend(sprintf('L*=%0.1f',L(checkers(1))/Lb(checkers(1))),...
        sprintf('L*=%0.1f',L(checkers(2))/Lb(checkers(2))),...
        sprintf('L*=%0.1f',L(checkers(3))/Lb(checkers(3))),...
        'FontSize',10);
    lgpos = get(lg,'Position');
    set(lg,'Position',[xp(i)+wd-lgpos(3) ...
        yp(i)+ht-lgpos(4) lgpos(3) lgpos(4)]);
end
for i = 1:8
    subplot(2,5,i);
    set(gca,'FontSize',12);
    grid on;
    ylim([0,1.5]);
    if i~=1 && i~=6
        set(gca,'YTickLabel',[]);
    end
    if i == 6
        ya = ylabel('PDF');
        set(ya,'Position',[-2.4 2 -1]);
    else
        ylabel('');
    end
    if i == 7
        xa = xlabel('v_y (mm/s)');
%         set(xa,'Position',[4.5e-3 -0.75 -1]);
    else
        xlabel('');
    end
end

cmat = parula;
subplot(2,5,10);
for i = [1,13,25, 26:28, 29,34,37, 38:43,...
        44,52,61, 62,69,78, 79,83,90]
    check = v(:,:,i);
    check = check(~isnan(check));
    bounds = minmax(check,[0.05 0.95]);
    check(check<bounds(1)|check>bounds(2)) = [];
    [bnnum,bnloc] = hist(check,numies);
    cval = round(L(i)/Lb(i) / max(L./Lb) * 256,0);
    cval = cmat(cval,:);
    plot(bnloc./max(bnloc),bnnum./max(bnnum),...
        '-','Color',cval,'LineWidth',1.5); hold on;
end
set(gca,'YTickLabel',[]);
xlabel('v_y/v_{y,max}'); ylabel('');
set(gca,'FontSize',12);
xlimies = [-1,1];
xlim(xlimies);
ylim([0,1.5]);
plot([-0.95 -0.8],1.3*[1,1],'-','LineWidth',2,...
    'Color',cmat(1,:));
tx = text(-0.75,1.3,sprintf('L*\\rightarrow0'));
set(tx,'fontsize',12);
plot([-0.95 -0.8],1.1*[1,1],'-','LineWidth',2,...
    'Color',cmat(end,:));
tx = text(-0.75,1.1,sprintf('L*\\rightarrow2.2'));
set(tx,'fontsize',12);
tx = text(0.15,1.32,'i) all (normalized)');
set(tx,'FontSize',14);

pause(1);
for i = 8:-1:1
    subplot(2,5,i);
    set(gca,'Position',[xp(i) yp(i) wd ht]);
end
subplot(2,5,10); set(gca,'Position',[x4 y2 wd*2+sp ht])


%% Cumulative density plots

Lb = (Kc./drho/g).^(2/3);
x1 = 0.10;
y1 = 0.76; y2 = 0.46; y3 = 0.13;
wd = 0.19; ht = 0.20; sp = 0.04;
x2 = x1+1*(wd+sp); x3 = x1+2*(wd+sp); x4 = x1+3*(wd+sp);
xp = [x1 x2 x3 x4 x1 x2 x3 x4 x1 x3];
yp = [y1 y1 y1 y1 y2 y2 y2 y2 y3 y3];
% colies = [0;0.6;0.85]*ones(1,3);
cmat = parula;

figure(60); clf;
set(gcf,'Units','normalized','OuterPosition',[0.1 0.1 0.7 0.9]);
for i = 1:8
    if i == 1
        checkers = [1,13,25];
        titlab = 'a) C1';
    elseif i == 2
        checkers = 26:28;
        titlab = 'b) C2';
    elseif i == 3
        checkers = [29,34,37];
        titlab = 'c) C3';
    elseif i == 4
        checkers = 38:40;
        titlab = 'd) C4';
    elseif i == 5
        checkers = 41:43;
        titlab = 'e) C5';
    elseif i == 6
        checkers = [44,52,61];
        titlab = 'f) S1';
    elseif i == 7
        checkers = [62,69,78];
        titlab = 'g) S2';
    elseif i == 8
        checkers = [79,83,90];
        titlab = 'h) S3';
    end
    subplot(3,4,i);
    xlimies = [0,0];
    for k = 1:3
        check = v(:,:,checkers(k))*1000;
        check = check(~isnan(check));
        bounds = minmax(check,[0.01,0.99]);
        check(check<bounds(1)|check>bounds(2)) = [];
        cval = round(L(checkers(k))/Lb(checkers(k)) / ...
            max(L./Lb) * 256,0);
        cval = cmat(cval,:);
        cp = cdfplot(check); hold on;
%         set(cp,'Color',colies(k,:),'LineWidth',2);
        set(cp,'Color',cval,'LineWidth',1.5);
        xl = [min(check),max(check)];
        xlimies = [min([xlimies(1),xl(1)]),...
            max([xlimies(2),xl(2)])];
%         pause(0.1);
    end
    xlim(xlimies);
    
    tx = text(xlimies(1)+range(xlimies)*0.70,0.65,titlab);
    set(tx,'FontSize',14);
    lg = legend(sprintf('L*=%0.1f',L(checkers(1))/Lb(checkers(1))),...
        sprintf('L*=%0.1f',L(checkers(2))/Lb(checkers(2))),...
        sprintf('L*=%0.1f',L(checkers(3))/Lb(checkers(3))),...
        'FontSize',10);
    lgpos = get(lg,'Position');
    set(lg,'Position',[xp(i)+wd-lgpos(3) ...
        yp(i) lgpos(3) lgpos(4)]);
end
pause(0.1);
for i = 1:8
    subplot(3,4,i);
    set(gca,'FontSize',14);
    title('');
    grid on;
    ylim([0,1]);
    set(gca,'YTick',0:0.5:1);
    if i~=1 && i~=5
        set(gca,'YTickLabel',[]);
    end
    if i == 5
        ya = ylabel('CDF');
        set(ya,'Position',[-90 1.3 -1]);
    else
        ylabel('');
    end
    if i == 7
        xa = xlabel('v (mm/s)');
        set(xa,'Position',[-3, -0.27 -1]);
    else
        xlabel('');
    end
end

cmat = parula;
subplot(3,4,[9,10]);
for i = [1,13,25, 26:28, 29,34,37, 38:43,...
        44,52,61, 62,69,78, 79,83,90]
    check = v(:,:,i);
    check = check(~isnan(check));
    bounds = minmax(check,[0.05 0.95]);
    check(check<bounds(1)|check>bounds(2)) = [];
    cval = round(L(i)/Lb(i) / max(L./Lb) * 256,0);
    cval = cmat(cval,:);
    cp = cdfplot(check./max(abs(check))); hold on;
    set(cp,'Color',cval,'LineWidth',1.5);
end
% set(gca,'YTickLabel',[]);
title('');
xlabel('v/v_{max}'); ylabel('CDF');
set(gca,'FontSize',14);
xlimies = [-1,1];
xlim(xlimies);
ylim([0,1]);
tx = text(xlimies(1)+range(xlimies)*0.04,...
    range(get(gca,'Ylim'))*0.88,'i) all');
set(tx,'FontSize',14);
plot([0.25 0.4],0.1*[1,1],'-','LineWidth',2,...
    'Color',cmat(1,:));
tx = text(0.5,0.1,sprintf('L*\\rightarrow0'));
set(tx,'fontsize',12);
plot([0.25 0.4],0.25*[1,1],'-','LineWidth',2,...
    'Color',cmat(end,:));
tx = text(0.5,0.25,sprintf('L*\\rightarrow2.2'));
set(tx,'fontsize',14);
set(gca,'Position',[0.15 yp(9) 0.3 ht]);

subplot(3,4,[11,12]);
xr = L./Lb;
% xr = PbB./Pf;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'o','markersize',8,...
    'MarkerFaceColor',[0,0,1],...
    'MarkerEdgeColor',[0,0,1],'LineStyle','none',...
    'LineWidth',2); hold on;
plot(NaN,NaN,'d','markersize',8,...
    'MarkerFaceColor',[0.4,0.7,0],...
    'MarkerEdgeColor',[0.4,0.7,0],'LineStyle','none',...
    'LineWidth',2);
plot(NaN,NaN,'d','markersize',8,...
    'MarkerFaceColor',[0.9,0.9,0],...
    'MarkerEdgeColor',[0.9,0.9,0],'LineStyle','none',...
    'LineWidth',2);
for i = 1:90
    plot(xr(i),cumMad(i)./cumMed(i),'Marker',symbols{i},...
        'Color',colors(i,:),'markerfacecolor',colorface(i,:),...
        'LineWidth',2); hold on;
end
% set(gca,'YAxisLocation','right');
xlabel('L*');
grid on;
ylabel('$\sigma/\overline{v}_{med}$',...
    'Interpreter','latex');
xlimies = [0,2.2];
xlim(xlimies);
ylim([0,4]);
tx = text(xlimies(1)+range(xlimies)*0.04,...
    range(get(gca,'Ylim'))*0.88,'j) stats');
set(tx,'FontSize',14);
set(gca,'FontSize',14);
leg = legend('C1','C2','C3','C4','C5');
set(gca,'Position',[0.57 yp(10) 0.3 ht]);
set(leg,'Position',[0.88 0.13 0.083 0.201]);

pause(1);
for i = [1,2,4,3,5,6,8,7]
    subplot(3,4,i);
    set(gca,'Position',[xp(i) yp(i) wd ht]);
    pause(0.001);
end



%%

xr = L./Lb;
yr = (Pf-PbB);
yr(yr<0) = 0;
yr = yr./PbB;
figure(105); clf;
for i = 1:length(xr)
    if i <= 43
        cval = colors(i,:);
        cvalf = cval;
    elseif i <= 78
        cval = [0.85 0.925 0.75];
        cvalf = 'none';
    else
        cval = [0.8 0.8 1];
        cvalf = 'none';
    end
    plot(xr(i),yr(i),'Marker',symbols{i},...
        'Color',cval,'MarkerFaceColor',cvalf,...
        'LineWidth',2); hold on;
end
set(gca,'XScale','log');
set(gca,'YScale','log');
xlabel('');
ylabel('');




%% Saving

if saveArg == 1
    save('Output.mat');
    
%     figure(10);
%     print(gcf,'Output\map_exchange','-djpeg','-r300');
%     figure(11);
%     print(gcf,'Output\map_vorticity','-djpeg','-r300');
%     figure(12);
%     print(gcf,'Output\map_velocity','-djpeg','-r300');
    
    figure(20);
    print(gcf,'Output\stat','-djpeg','-r300');
    
    figure(30);
    print(gcf,'Output\mod_v1','-djpeg','-r300');
%     figure(31);
%     print(gcf,'Output\mod_v2','-djpeg','-r300');
%     figure(32);
%     print(gcf,'Output\mod_v3','-djpeg','-r300');
    
    figure(40);
    print(gcf,'Output\mod_u','-djpeg','-r300');
    
    figure(50);
    print(gcf,'Output\mod_vort','-djpeg','-r300');
    
    figure(60);
    print(gcf,'Output\CDF','-djpeg','-r300');
%     figure(61);
%     print(gcf,'Output\PDF','-djpeg','-r300');
end


%% Function

function mod = testeq(fa,fb,fc,fx)
    mod = exp(fa*-fx.^fb)+fc;
end
function mod = testeqvi(fa,fb,fc,fd,fx,fy)
    mod = fa*exp(-fb*fx.^fc)./fy.^fd;
end

